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Abstract 

The quantum modes of a nonlinear Klein Gordon lattice have been computed numerically [L. 
Proville, Phys. Rev. B 71, 104306 (2005)]. The on-site nonlinearity has been found to lead 
to phonon bound states. In the present paper, we compute numerically the dynamical structure 
factor so as to simulate the coherent scattering cross-section at low temperature. The inelastic 
contribution is studied as a function of the on-site anharmonicity. Interestingly, our numerical 
method is not limited to the weak anharmonicity and permits to study thoroughly the spectra of 
nonlinear phonons. 

PACS numbers: 63.20.Ry, 61.10.Dp, 61.12.Bt, 61.14.Dc 
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I. INTRODUCTION 



The metal hydrides are typical compounds whose the technological importance^ aroused 
many spectroscopic inquiries. The spectroscopy of hydrogen modes allows to work out the 
occupied interstitial gaps, e.g. octahedral or tetrahedral and thus the phase structure of 
the alloy. Further, the determination of the anharmonicity of hydrogen sites^»^ permits 
to evaluate the on-site potential landscape'^ and so the hydrogen ability to migrate or to 
dimerise. The inelastic neutron scattering (INS) revealed the on-site anharmonicity in the 
model hydrides as NbH, TaH and PdH (for a review see Ref. m). In the compounds TiH 
and ZrH, that anharmonicity was found^i^ to lead to the optical phonon bound states, i.e., 
some well-known excitations in molecular crystals^'^'^'^iii (also quoted as biphonon and 
bivibron) where the anharmonicity of internal covalent bonds overpasses the interaction 
between molecules. As in the metal hydrides, the proton dynamics has been studied by 
INS in molecular crystals as polyglycine^^ and 4-methyl pyridine^^. In the latter, the bound 
states of the methyl group rotational modes proved to last several days^*^ (see also Refs. 

The trapping of energy, whether it is in proton vibrations or in the intrinsic modes of 
molecules is a consequence of the emergence of breather, a theoretical paradigm that has been 
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intensive research in different contexts (for instance see Refs. 
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29l ). The distinctive property of breather is to gather the spatial localization 



and time periodicity. For the case of a translational-invariant quantum lattice, the phonon 
bound states can be thought as the siblings of breathers^'^'^'^, as phonon bound states 
and breathers both stem from the lattice anharmonicity. The formers appeared, though, 



earlier in literature (see Refs. 



25 



321 and references therein). The link between breather 



and phonon bound states has been studied in different work9i^i2ii2^i2ii2i. Consequently, the 
experiments that enhanced the existence of phonon bound states (for instance see Refs. 
^I^ IO) are equally some concrete evidences of the emergence of quantum breathers. In 
addition, several experiments^'^'^ were recently dedicated to breather and its important 
role in the energy storage and transport. As an overview, one may retain that the studies 
on breather aim at improving the standard harmonic treatment of lattice dynamics which 
flaws are well-known in solids—. 

The scattering of neutrons^E^ X-rays^ and electrons^^ provided numerous insights into 
condensed matter and molecular physics. In particular, the inelastic scattering allows to 
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probe the dynamics of crystals and molecules. The dynamical structure factor, S{q, uj) which 
is proportional to the scattering cross-section can be calculated by solving the Schrodinger 
equation for the lattice modes (within Born approximation), so this quantity is essential for 
bridging theory to experiment. However the lattice eigen-modes may be figured out exactly 
only in the harmonic lattice model (see for instance Ref. |39|]42| ) and thus the anharmonic 



contributions to the energy are neglected, although they stem from the atomic interaction 
potential. In order to examin the effect of anharmonicity on S{q,uj), some approaches have 
been attempted in different quantum lattice models, e.g., the Hubbard model for bosons^i^ 
and the Klein Gordon model with a weak on-site anharmonicity^. The latter is often quoted 
as nonlinear Klein Gordon lattice (KG) and, noteworthy, it may account for the intrinsic 
anharmonicity of molecular crystals^i^^^iii, molecule bonds^i^ or metals where interstitial 
gaps are filled with light particles, as metal hydrides^'^'^'^. It is the reason why we proposed 
a numerical method^ that is tractable for different type of nonlinearity, to compute the non- 
linear quantum modes. In the present paper, our previous developments are used to evaluate 
the coherent Dynamical Structure Factor (DSF) of the KG lattice, at low temperature. The 
contribution to the DSF of the nonlinear quantum modes, known as either the phonon bound 
states^'^'^'^i'^ or else the quantum breathera^^'^^'^, is the central purpose of our work. As 
the DSF standard derivation is obtained in the harmonic approximation^i^, treating the 
anharmonicity as a perturbation, we propose a different scheme where this approximation 
is not required. For instance, the Bloch identity, which is a key for the conventional calcula- 
tion, is not invoked in our theory. We simply introduce a Taylor expansion of the DSF with 
respect to the atomic displacements, computing numerically the coefficients of the series. 
Our method is tested by comparison to the standard analytical calculation on the purely 
harmonic lattice. Our approach allows us to deal with different strength of the anharmonic- 
ity. When the on-site nonlinearity dominates the intersite coupling, the nonlinear quantum 
modes lead to some anharmonic resonances well separated from the multi-phonon continua, 
in the inelastic spectrum. These anharmonic resonances have same order as for fundamental 
phonons whereas the unbound multi-phonons yield a DSF at least one order below. On the 
other hand, their dispersion is found to decrease dramatically with the value of the energy 
transfer, which differentiates them from the broad phonon resonance. By way of contrast, 
when the intersite coupling is larger than the on-site nonlinearity, the anharmonicity of the 
spectrum is weak compared to the phonon dispersion. In such a quasi-harmonic lattice. 
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the dynamical response of the biphonon (a two phonon bound state) may yet exhibit a 
significant magnitude, still much larger than the two-phonon, even though the biphonon 
resonance occurs at an energy transfer that approaches the unbound phonons band. Actu- 
ally, the biphonon signature is found to dominate the two-phonon DSF provided that the 
lattice nonlinearity is not strictly zero. Alongside our numerics, a perturbation theory is 
developed and proves accurate for strongly nonlinear lattices. We also show how the model 
parameters may be adjusted so as to simulate the inelastic scattering of a material and thus 
to work out the energy landscape of the inner particles as well as their interactions. Finally, 
the present study will serve as a basis for future simulations of the KG incoherent scattering 
cross-section. 

After a brief introduction to the nonlinear KG lattice, the DSF is derived in Sec. HIl 
in Sec. IIIII Results and comments are detailed in the same section while the Sec. IIVI 
summarizes them and draws some perspectives. 

II. COHERENT DYNAMICAL STRUCTURE FACTOR OF THE NONLINEAR 
KG LATTICE 

We assume that some light particles form a regular network, whether it is inside a molecule 
or a crystalline solid. Were the particle dynamics independent, the Hamiltonian would have 
read: 



where xi and pi are displacement and conjugate momentum (i.e., [a;/,pz] = hi) of a particle 
at site /. For sake of simplicity, we choose to work on a one-dimensional lattice with a 
single direction for the motion of atoms, denoted by the unit vector u, which is named 
as polarization in the foUowings. The single dimension proved relevant in several concrete 
studies^'^'i^i^. The local potential V{xi) may then be expanded as a fourth order series: 
V{xi) = a2xf + a^xf + a^x^. Some higher order terms could be added with no difficulties for 
our numerics, so as to simulate eventually a specific shape of V . In truth, the interaction 
between nearest neighbors, that may be direct or mediated by the heavy atoms that compose 
the skeleton of a material, involves a displacement coupling that is modeled by a quadratic 




(1) 
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term. So the Hamiltonian of the particle ensemble now reads: 

H = Ho-c i^i-^jf, (2) 
i,j=<i> 

where c is the coupling parameter and the symbol (< Z >) describes the first neighbors of 
the site labeled by /. Introducing the dimensionless operators Pi and Xi, the Hamiltonian 
can be reformulated as follows: 

I j=<i> 
where the fundamental frequ ency Q, as well as the dimensionless coefficients A^, and 



C have been defined in Ref. 
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The total number of sites is denoted by A^, the lattice 
parameter by Oq and the orientation of the chain is defined by a unit vector v. We introduce 
ip, the eigenstates of H, and particularly the groundstate ipcs, the phonons ipphi^) and the 
biphonons ipbiiX)- These states are computed numerically^ as well as their eigenvalues Eos, 
Eph(k) and Ebi(k), respectively. The wave momentum, denoted by k, verifies k = 
where k is an integer. The groundstate wave momentum ko takes its value in the reciprocal 
lattice, ko = where is an integer. Apart when it will be noted, ko is fixed to 

zero. The lattice is assumed to be maintained at very low temperature, i.e., kT < Eph, 
so the scattering induces some excitations solely from ipcs- The transition probability is 
proportional to the DSF— : 

5(q,^^(k)) = ^1 E < fele^^^-'-^ll^(k) > r (4) 

j 

where ci;,/,(k) = (-Si/)(k) — E^gs)/^^ ^ the scattering vector and rj represents the atomic 
position at site j. That position can be expressed as rj = jaov+CXjU, where C = ^Jfij (mil) 
is the length scale of vibrations, fixed to a reasonable value of 3% of oq. The polarization 
vector u is the principal axis of atomic displacement around the equilibrium position (jaov). 
We now spell out our method for computing the DSF. As, strictly speaking the Bloch 
identity does not hold for a nonlinear lattic^^, we propose a derivation which differs from 
the conventional one. Since C is much smaller than oq, the exponential function in Eq. (jl]) 
can be expanded as a Taylor series with respect to the weakest of its arguments, the term 
proportional to C: 

^(q,^.«) = ^1 E^^^"^^"" X E < tel^JI^(k) > r (5) 



This development holds provided that the order of the expansion is large enough. The series 
convergence has been tested by increasing the order up to reach the desired precision. The 
agreement to the standard analytical calculation on the purely harmonic lattice (see Sec. 
imp is a necessary condition of validity that has been fulfilled too. In the left hand side of 
Eq. the bracket < ipGsl^jli'iX) > can be replaced by < ipcsl^ol'^O^) > multiplied by a 
phase factor e"**^''' "^]"''''^ because of the translational invariance. The sum over the subscript 
j in Eq. ([5]) gives zero for all q, aside from the wave vectors that match the momentum 
conservation [(q — k).v] = 0. This condition, as well as the conservation of the energy are 
supposed to be satisfied for a given state ip{k). The calculation of ^(q, cu^) can then be 
reduced to the determination of the bracket D{ip(k),p) =< ipcsl^ol'^i^) since one may 
check that: 

5(q,o;^(k)) = iV^[(q-ic).v]|^ ^'^'^-'|^^^' /^(^(k),p)r (6) 

p 

We study the variation of the DSF along a single direction which the angles with respect to 
u and V are and a^, respectively. Then S'(q, cj^) depends only on the magnitude of the 
momentum transfer |q| along that direction and the conservation law fixes |q| cos (a^) = [k.v]. 
This condition can be achieved whether \a^\ 7^ 7r/2 which may be reasonably assumed for a 
low-angle scattering. In a same manner, when = tt/2 the inelastic part of the structure 
factor is identically null because [q.u] = (excepted for the elastic DSF which then equals 
N). Consequently, we assume that \ay \ = \au\ 7^ 7r/2 which keeps the physics of our problem, 
avoiding the unimportant cases. We denote by q the scalar product [q.u] for a vector q that 
matches the conservation law, so that we have also q = [k.v]. Then, the DSF reads: 

S{q,u;^^,^) = N\Y,^^D{^|J{q),p)\\ (7) 

Since the optical modes are characterized by a weak value of C in Eq. ([3]), the first step of 
our treatment is concerned with a perturbation theory with respect to the intersite coupling. 
At C = 0, namely the anti-continuum (e.g. , uncoupled, atomic or molecular) limiti^, the 
phonon states and the phonon bound states can all be written as some Bloch waves^: 

where (paj is a on-site wave function that depends only on Xj. Actually, (paj is the ath 
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eigenstate of the on-site Hamiltonian: 

hi = ^ + A^Xf + A,Xf + A,Xf. (9) 

To a first order in C, the lattice groundstate iIjgs is simply given by the product Bq = Ili4>o,i 
while Eq. (jHl) gives a phonon for a = 1 and a biphonon for a = 2. Let us denote by Va,n, 
the projection of the state (pa onto the nth on-site harmonic oscillator eigen-state \n > and 
develop the bracket G{a,0,p) =< (f)o\XP\(f)a >■ 

G(a,0,p)= Yl Vo,n.Va,n<m\XP\n> (10) 

n>0,m>0 

where the subscript has been dropped for ease. We point out that for (7 = 0, we have 
D{Ba(k),p) = G{a,0,p)/VN and D{Bo,p) = G{0,0,p). One more thing has to be done 
before reaching our goal, is to compute < m\XP\n >. To that purpose, the Bose-Einstein 
operators [i.e. a"*" = -\/2(X — iP) and a = V2{X + iP)] are used to expand X^. We wrote a 
fortran program which realizes the expansion, respecting the commutation rule, [a, a~^] = 1. 
For example, the output for p = 10 is: 

= l[a+io + aio + io(a+a9 + a+9a) 

+ 45(a+^a^ + a+^a^ + + a^) 

+ 120(a+^a^ + a+V) + 360(a+a^ + a+^a) 

+ 210(a+V + a+^a^) + 1260(0+^0^ + a+^a^) 

+ 630(a+^ + a^) + 252a+^a^ + 2520(a+^a^ + a+^a^) 

+ 3780(a+^a + a+a^) + 3150(a+^ + + a+^a^) 

+ 4725(2a+2a4 ^ 20+^0^ ^ 4^+2^2 ^ ^+2 + ^2 ^ 2a+a) 

+ 12600(a+a^ + a+3a + a+V) +945]. (11) 

The writing of X^ would take more than one full page for values of p larger than 30. Our 
program allows us to compute the bracket < mlX^'Ira > up to the power p = 60, for different 
integers m and n. The association of this program with the numerical diagonalization of hi 
(Eq. (^) which fixes the coefficients V^^„ in Eq. flTIl — . permits to compute the coefficients 
G{a,0,p) that can be tabulated for different model parameters. Finally, to a first order in 
G, one obtains: 

p 
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for the inelastic scattering where the momentum conservation imposes q = 2'n'k/{Nao), 
whereas the elastic response is given by: 

5(g,0)=iV|5^^^G(0,0,p)p (13) 

where q = 27rA;o/cto [k and range over integers]. Formally, these results are given for a 
one-dimensional lattice but they can be extended to higher dimensional lattice by summing 
over the coordinates and polarizations. To a first order in C, the frequency ujBa{q) can also 
be evaluated by: 

ujB^iq) = ^{-fa - 70 - 2C X G{a, 0, 1)^ cos(gao)) (14) 

where the coefficient 7q, is the eigen-energy of hi (Eq. ([9])) associated to $„. For the perfect 
harmonic lattice, the first order in C in the standard calculation^^i^ of S{q,uj) is equivalent 
to Eqs. ([IID and ([T3D. 

When C is large, the Hamiltonian in Eq. ([3]) must be diagonalized numerically after 
expanding in a suitable basis. We worked with a Bloch wave basis^ given by: 

i?[n,«,](k) = X nz<^„,,+, (15) 

where A^n^^aj] ensures the normalization. Each eigenstate, characterized by a wave vector k, 
can be written as a linear combination of those Bloch waves: 

V.(k) = 5^W^^,n,„,i?[n,„,](k), (16) 

where the subscript HjCtj identifies a single Bloch wave in Eq. f|T5|) . The numerical diagonal- 
ization has been carried out for different lattice sizes with no noticeable discrepancy in the 
eigenspectrum in increasing N. In Fig. [H the same energy cutoff on the Bloch wave basis 
has been fixed for both = 23 and = 33. Apart on their wave vector, the eigenvalues are 
found to be independent of A^. In Fig. [1] (a), the optical phonon branch is plotted whereas, 
in Fig. [1] (b), the energy region of the first overtone is plotted for same parameters as in 
(a). One notes clearly the biphonon branch and the two-phonon band, thoroughly described 
elsewhere^'i^i^'^'^. Interestingly, for small enough nonlinear parameters, the anharmonicity 
of the lattice is negligible compared to the phonon branch width and the biphonon branch 
disappears (see Fig. [1] (c)). The bracket D{ip,p) can now be written as follows: 

D(^(k),p) = Yl ^;Gs,n,/3,^^.n,«, < 5[n,ft.](0)|XonS[n,„,](k) > (17) 
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where the last term in the right hand side can be detailed further: 

xUi-^j6{ai+ni,f3i+n2)G{ai+n2,<^i+ni,p)- (18) 

This equation concludes our computation task on the KG lattice dynamical response to 
the scattering of a external beam of light or particles. Our numerical treatment, including 
the diagonalization of H and the computations of Eqs. U\ [13 and [18] takes few hours on 
a conventional desktop PC. The convergence of the series in Eq. has been tested by 
comparing our results for different development orders, e.g. 30, 40, 50 and 60. As expected, 
the smaller is the scattering vector, the better is the precision. For a momentum transfer 
that ranges over 30 lattice Brillouin zones, the difference between the DSF computed with 
series orders 50 and 60 is less than 1 %. 



III. RESULTS AND DISCUSSIONS 

The Fig. [2] shows our typical result that is a 3D plot of the inelastic S{q,uj) for a 
nonlinear lattice which model parameters are those of Figs. [1] (a) and (b). Three noteworthy 
resonances emerge with same order of magnitude. These three resonances correspond to the 
eigen-energies of the nonlinear phonon states, i.e., the phonons, biphonons and triphonons. 
The binding energy of the biphonon, evaluated at the center of the Brillouin zone (BZ) is 
around a few percent of the fundamental phonon excitation, as measured in different metal 
hydrides-'*'. Several model parameters have been tested and lead to qualitatively similar 
results. For instance, besides the different values of the energy transfer, the biphonon and 
triphonon resonances occur in a lattice with a quadratic-quartic on-site potential, i.e., = 
in Eq. [3l In the biphonon and triphonon resonances, one may recognize some peaks that will 
be examined further. In order to analyze our results more quantitatively it is, though, easier 
to work with a 2D plot that represents the S{q,u!) profile versus the projection q = [q.u] 
of the scattering vector q. In order to verify the accuracy of our computations, presented 
in Sec. [TTl we compare our numerical results to the exact analytical ones that are achieved 
in a purely harmonic lattice. In that case, the profile of the Debye- Waller factor, i.e., 
S{q, 0)/N is plotted in Fig. [3] and shows a convincing agreement since our numerical results 
scarcely differ from the analytical ones. It ensures us that our series expansion of S{q,u) 
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has converged. Only the non-zero values of the Debye- Waller have been retained in the 
plot, i.e., q = 27rA;o/ao. A similar accord is obtained for the inelastic part of S{q,u!) (see 
further in the same section). The calculation made in a standard harmonic approximation^ 
gives a Gaussian dependency of the Debye- Waller factor S{q,0)/N = exp{—2W{q)) where 
2W{q) = {qC)'^\ < 4'Gs\^j\4'GS > I- In this expression, the bracket is simply given by 
V(2iV)Exl/^p/^(^) where UphiK) = ^1 + 2Ccos{K.ao) and K ranges in the first BZ. 
When the nonlinear parameters are no longer negligible in the Eq. ([3]), we note, in Fig. [3] 
that the Debye- Waller of a nonlinear lattice differs from the harmonic one. The Gaussian 
form is yet roughly conserved so there is no qualitative changes involved by the on-site 
anharmonicity. Our perturbation theory proves sufficient under the condition that the inter- 
molecular coupling is not too large. In agreement with our results, the earlier study of 
B.V. Thompson^ concluded that a cubic-anharmonic term involved a negligible correction 
upon the Debye- Waller factor. The calculation of Thompson is, in fact, a second order 
perturbation in the cubic term, which is thus assumed to be weak enough for the perturbation 
theory to be valid. Not surprizingly, it leads to a correction that is also weak. However, as 
we shall see, the most significant contribution of the nonlinearity to the DSF spectra proves 
to be inelastic and is due to the phonon bound states. 

In a same manner as for the elastic scattering, our numerical computation of the in- 
elastic DSF is compared to the exact calculation^ in the case of a harmonic lattice. The 
contribution of phonon states is then given by: 



with same notations as previously used in the same section. The obtained agreement, 
demonstrated in Fig. HI confirms the validity of our theoretical approach and consequently 
the convergence of our series development. The DSF profile appears as being continuous 
because the lattice size is large enough to blur the discreteness of the Fourier space at the 
scale of Fig. HI which ranges over 60 BZ. Apart from the S{q,uj) ripple, the perturbation 
theory captures quite well the main variations of the envelop of S{q, ujph)- This ripple stems 
from the dispersion of Uph and so, its amplitude increases with C to become some peaks for 
the acoustic phonons in a harmonic lattice^. When the sign of C is inverted in our model 
(C becomes negative), the local minima of S{q,ujph) are shifted at the edges of the lattice 
Brillouin zones instead of being in the middle when C > 0. In Fig. |U the DSF resonance 




(19) 
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involved by the two-phonon states is also plotted. It is usually termed multiple-scattering 
and is one order of magnitude smaller than the one-phonon process. If the anharmonicity 
of V can not be neglected, the profile of the one-phonon resonance in S{q, u) shows no 
qualitative changes compared to the harmonic lattice (see Figs. IHandE] (a)). However, the 
response of the biphonon, which profile is denoted by S{q,ujhi), appears clearly at higher 
energy transfer than one-phonon (see Figs. [2] and [5] (a)). The resonance associated to the 
biphonon is one order of magnitude larger than for the unbound phonon states and has same 
order as S{q, ujph)- A gap is opened between the biphonon branch and the two-phonon band 
(Fig. [T] (b)) so the biphonon resonance occurs for an energy transfer which does not match 
the harmonics of the fundamental phonons. As mentioned above, that energy gap is also 
called the binding energy of the biphonon^. It is also worth noting that the maximum of the 
S{q, uJbi) envelop is reached for a scattering vector that differs from the envelop maximum 
of S{q,u!ph)- This maximum occurs, indeed, at larger q for the biphonon. The S{q,uJhi) 
envelop can be evaluated within the perturbation theory (see Figs. [5] (a)) which provides a 
rather satisfactory approximation, although the S{q,uJbi) ripple has a larger amplitude than 
for phonons. To interpret our results, it is needed to dwell on the ripple of the biphonon 
contribution to the DSF. As for the phonon, this ripple is related to the dispersion of the 
biphonon branch. In the case where the binding energy of biphonon is not very large, say 
no more than a few percent of the phonon energy, the different bottoms of the S{q,uJhi) 
ripple occur at center of the successive BZ while the tops are reached at the BZ edges. 
According to several tests, it is a systematic feature which, in contrast to phonon, does 
not depend on the sign of C. Moreover, as shown in the insets of Figs. [5] (a) and (b), 
the minima of S{q, uJu) and the maxima of the two-phonon resonance correspond one to 
one. That may be explained as follows. In our perturbation theory the biphonon states are 
approximated by the Bloch waves that bear a single on-site excitation, a = 2 in Eq. ([8]). 
In case C 7^ 0, since the biphonon binding energy is not very large, these Bloch waves are 
hybridized with some other Bloch waves, given by Eq. (|T5i) . and particularly those bearing 
two on-site excitations a = 1, at different sites. The degeneracy-lifting of the latest states 
yields the unbound two-phonon band. The ripple of the DSF can not be analyzed through 
our first order perturbation theory (see Figs. [5] (a) and (b)), which hints that it is due to 
the Bloch waves hybridization, involved by the intersite coupling. The discrepancy of the 
perturbation theory Eq. (|T2l) increases as the biphonon gap decreases (compare Figs. [1] (b) 



11 



and (c) to Figs. [5] (a) and (b)) because of the hybridization step-up. The gap between the 
biphonon branch and the unbound two-phonon band is smaller at the center of each BZ (see 

n 

Fig. [T] (b) and Ref. 311) because here the width of the two-phonon band is maximum. The 
smaller the energy gap is the larger the hybridization is, so as a result, one finds a larger 
contribution of the two on-site excitation Bloch waves into the biphonon eigenstate at the 
center of each BZ. The Bloch waves with multiple excitations at distinct sites yield a zero 
response to scattering, S{q,uj) = 0. One deduces that the biphonon response S{q,Ubi) is 
minimum when the contribution of the two on-site excitation Bloch waves is maximum, i.e., 
at the center of each BZ. On the other hand, the unbound two phonons resonance comes to 
a maximum at the center of each BZ because of the contribution of the Bloch waves with 
a single on-site excitation a = 2 to the two-phonon eigenstates. Conversely, at the edge of 
the lattice Brillouin zones, the Bloch wave hybridization is minimum (because the energy 
gap is maximum) so that the biphonon response is maximum and may even reach the value 
computed within the perturbation theory (Fig. [5] (b)). Another interesting point that is 
made clear within the previous discussion is why the S{q,uJbi) ripple is not shifted when 
the sign of C is changed, unlike S{q,u!ph). Indeed, reversing the C sign leaves unchanged 
the two-phonon band shape, i.e., the band width is still maximum at the BZ center. So it 
is for the biphonon branch which the dispersion is mainly due to the contribution of the 
two on-site excitation Bloch waves^i (the dispersion computed within Eq. f|T^ for a = 2 is 
much smaller than the effective biphonon band width). So the gap between the biphonon 
and two-phonon bands is equally unchanged under the inversion of the coupling sign, i.e., 
the energy gap is still maximum at the BZ edges and minimum at the center. According 
to our analysis of the S{q, uju) ripple, it makes clear why the maxima and minima of the 
biphonon DSF are independent of the sign of C. Our argumentation holds provided that the 
biphonon binding energy remains inferior to a few percent of the fundamental excitations. 

When the anharmonicity of V is weak compared to the phonon band width, the biphonon 
gap closes and a pseudogap formal at the edge of the lattice Brillouin zones. Then the 
biphonon branch hardly appears nearby the two-phonon band (see Fig. [D (c)). The lattice 
is said quasi-harmonic since the spectrum anharmonicity vanishes. There is, indeed, no 
anharmonic resonances, neither in the density of state nor, accordingly, in the DSF. That 
may happen for higher order phonon bound states as the hardening of the quartic on-site 
term in V (Eq. ([3])) could compensate the softening of the cubic potential in a certain range 
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of energy. Although that seemingly perfect harmonicity, the S{q, u) ripple may yet betray 
the nonlinearity of V. In truth, at the BZ center, the DSF associated to the biphonon, 
referred to as S{q,ujf,i), has same order as the DSF of two-phonon states, which magnitude 
is one order below the phonon resonance (see Fig. El^b)). In contrast, the profile ^(g, cj^j) 
reaches its maxima where the pseudogap opens (see Fig. [U (c)), i.e., at the BZ boundaries 
and the ripple alongside each BZ can be seen as a series of peaks. In Fig. [5] (b), the largest 
peak remains much larger than the two-phonon response and it has same order as the 
maximum of the phonon DSF profile, S{q, ojph). The maxima of the biphonon DSF decrease 
as the nonlinear parameters tends to zero. However, before to reach that point, we have seen 
that the anharmonicity of the energy spectrum vanishes as in Fig. [T] (c). Consequently, we 
propose to dub the S{q, Ubi) peaks by nonlinear Bragg peaks, since they gather two features 
that are opposite to the usual Bragg scattering: first, they are inelastic peaks and second, 
they satisfy the Bragg refiection condition but shifted of a half reciprocal lattice vector 
(i.e., vr/ao for the one-dimensional chain). According to Ref. I47I, the pseudogap may occur 
in lattices with a higher dimension so that similar results as those presented here can be 
expected in different lattices. In spite of the relative simplicity of our model compared to the 
case encountered in practice, it indicates that the biphonon may still emerge and contribute 
significantly to the scattering, even though the spectrum seems perfectly harmonic. In Figs. 
[T](c) and[5](b), we have chosen intentionally a set of parameters that yields a eigen-spectrum 
which is almost harmonic. A pseudogap hardly opens between the biphonon branch and 
the two-phonon band. The parameter could be doubled [as in Fig. 2 (c) in RefQ 
without opening a substantial gap which would permit to identify the anharmonicity in the 
energy density of state. In fact, the larger is the phonon dispersion the larger is the range of 
nonlinear parameters that leads to a quasi-harmonicity. In Fig. [6] (a) the profile of S{q,uj) 
is plotted for different values of the energy transfer u, e.g., Uph for phonon, for biphonon 
and ujtri for triphonon. The triphonon resonance is found to reach a significant magnitude 
that is comparable to phonon and biphonon. Such a high order phonon bound state has 
been identified in different materials, e.g., TiH^ and ZrH^ and HCl^S. It is worth noting the 
triphonon spectrum ripple, which originates from the hybridization between the bound and 
unbound phonon states, as for the biphonon spectrum. This triphonon ripple is shifted of 
vr/ao compared to the biphonon DSF profile. 

In what follows, we spell out how our theory could be linked to some concrete cases. 
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To simulate a scattering experiment, we must account for the resolution of the technics. If 
one assumes that the energy resolution is not sufficient to differentiate the inelastic lattice 
resonances (as it may be the case, for instance, in the inelastic X-rays scattering), we have 
to integrate our simulation with respect to the energy transfer uj. The elastic contribution is 
skipped from that integral so as to distinguish the inelastic part of the DSF. We computed 
the corresponding integral and reported the result as a dashed line in Fig. [6] (a). The 
amplitude and the position of the maximum of the DSF integral are related to the nonlinear 
parameters. In a strictly harmonic lattice, that integral scarcely differs from the phonon 
response, whereas in Fig. [6] (a) we note a clear difference between the dashed line and the 
profile S{q,ujph). The envelop of the DSF integral has a form that might be fitted by a 
superposition of two Gaussian functions, centered at different momenta. It is a standard 
treatment in the interpretation of the experimental S{q, uj) profile (see for instance Fig. 8 



in Ref, 



13l ). If the nonlinear Bragg peaks are large (i.e., when the anharmonicity is weak 



compared to the phonon dispersion but that is yet not negligible) they may emerge in the 
integral as shown for large q in Fig. [6] (a). In a strongly nonlinear lattice, a large gap 
separates the biphonon branch from the two-phonon energy band so the nonlinear Bragg 
peaks shrinks to a ripple. Then it becomes rather difficult to depict the contribution of 
the nonlinear states in the DSF integral, aside from the shift of the envelop to larger q. 
On the other hand, if the energy resolution is sufficient to resolve the energy dispersion of 
the phonon branch but the accuracy over the scattering vector is not to separate the BZ 
(as it may be the case, for instance, in the neutron scattering of a powder specimen, since 
u and V are randomly distributed), then our simulation must account for the contribution 
of the distinct BZ. As in a one-dimensional lattice the only symmetry is the inversion, we 
have to sum the factor 2S{q + ko,u{q)) over the reciprocal vector ko, where q and uj{q) are 
both assumed to match the momentum and the energy conservation, respectively. Since 
S{q + ko, uj{q)) decreases exponentially for a large enough value of (g + ko) (see Fig. [6] (a)), 
the sum is ensured to be finite. The result of our treatment is shown in Fig. [6] (b). We 
note the three resonances due to the nonlinear states, i.e., the biphonon, the triphonon and 
the quadriphonon that are separated from the broad resonances of either the phonon or the 
unbound phonons. Here, we must enhance that whether the inelastic scattering is coherent 
or incoherent, the lattice resonances occur at same energy transfer. Thus, from that point 
of view it is relevant to compare our simulation to the experimental spectrum of incoherent 
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INS in metal hydrides. Actually, the model parameters in Fig. [6] (b) has been adjusted so 
as to obtain the same energy resonances for the phonon, biphonon and triphonon as in the 
spectrum shown in Fig. 2 of Ref. |49|. The energy unit is HQ ^ 107meV, the dimensionless 
coupling is C = 0.06 and the nonlinear parameters are = 0.17 and ^4 = 0.0231. In that 
manner, the on-site nonlinearity and thus the energy landscape of the light particles, here 
the deuterium, may be worked out to a better accuracy than in a purely quadratic model 
(see the inset in Fig. [6](b)). The same simulation could have been achieved in other alloys as 
TiH and ZrH^i^. In their earlier studies, A.I. Kolesnikov and co-workers^i^i^ recognized the 
phonon bound states resonances in the INS spectra of TiH-D and ZrH-D. To analyse their 
experimental spectra, the authors compared their data to a simulation made with a model 
proposed initially by V.M. Agranovicb^, i.e., a one-dimensional boson Hubbard lattice. In 
such a model it is, though, difficult to figure out the potential landscape of the interstitial 
gaps. The single dimension of the model was shown to be reasonable because of the strong 
anisotropy of the H-H interaction along the metal c-axis. Since we are studying the coherent 
DSF, a accurate simulation of the spectra obtained by Kolesnikov et al. for the incoherent 
INS is out of purpose. We shall attempt such a exercise in a future work, devoted to the 
computation of the incoherent DSF. The amplitudes of various resonances as well as their 
widths might then be inferred. In Fig. [6] (b), our simulation has been realized with a typical 
energy resolution function, i.e., a triangle function centered at uj with a width of 0.02ti;, which 
corresponds to the crystal analyser TFXA, ISIS Rutherford Appleton Laboratory, described 
in Ref. [gI. The dynamical trajectory of the TFXA machine has not been simulated in Fig. 
6 (b). It is noteworthy that the amplitude of the nonlinear resonances has same order as 
for phonons. Similar calculations with different model parameters even showed that, for 
a stronger anharmonicity the nonlinear resonances may be even larger in magnitude. The 
reason for such a behavior is that the density of state enhances as the eigenstates band 
width decreases. This effect is sharp around the narrow branches as those of phonon bound 
states^ (see Fig. [1] (b)) which explains why the biphonon resonance may dominate the 
phonon one. The increase in the density of state may be recognized too in the single phonon 
response because the phonon branch is flat at the edge and at the center of the lattice BZ 
(see Fig. [T] (a)). Thus one sees two side-bands in the phonon resonance, at the upper and 
lower boundary in energy. These side-bands may be identified as phonon wings. According 
to our reading of the above cited works, the experimental INS spectra of the metal hydrides 
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do not exhibit that property but, once again, one must left that point aside, until we study 
the incoherent DSF. On the other hand, it is worth noting that the form of the lattice 
response around the first and the second overtone energy regions, in Fig. E] (b), is similar 
either to the infrared adsorption spectra^ in crystals as CO2, N2O and OCS or to the Raman 
spectrum of H2 solid^^. Indeed the emergence of a sharp peak, associated to a biphonon (or 
to a bivibron^^) occurs near a small hill-like resonance that is due to the unbound phonon 
states. In our simulation, a gradual variation of the model parameters so as to decrease the 
anharmonicity shows the same behavior as the pressure-induced bound-unbound transition, 
at 25 GPa in II2 solid^, or at 34 GPa in D2 solict^. Around that transition, the biphonon (or 
bivibron) peak broadens and decreases in magnitude. In our model, the pressure variation 
can be simulated by a change of the coupling parameter C due to the fact that neighbouring 
molecules are moved closer together because of the external pressure. 

Although our work is concerned with a one- dimensional lattice, the results enhanced in 
the present study should hold in higher dimension where the DSF depends on the lattice 
symmetries and polarizations. The extension to higher dimension could be achieved in our 
perturbation theory with no particular difficulty. This approach proves adequate under the 
condition that the intersite coupling is very weak, as expected in the metal hydrides as NbH 
and TaH, from the study of J. Eckert et al?. Our numerics, more accurate, could also be 
extended to higher dimension but that would require either to use a powerful computer or 
else a dramatic restriction on the site number, which could yet be relevant for some molecules 
as benzen^. 

IV. CONCLUSION 

The dynamical structure factor (DSF) of the nonlinear Klein Gordon chain (KG) has been 
calculated for different strength of the on-site anharmonicity. This has been possible thanks 
to our numerics that permit to diagonalized accurately the non-quadratic Hamiltonian^. 
The DSF has been expanded as a Taylor series of the atomic displacements, which avoids 
the use of the conventional harmonic approximations. We found that the on-site nonlinearity 
leads to phonon bound states and consequently to some anharmonic resonances in the DSF 
spectrum, which somehow confirms other works on different quantum lattices^'^'^. The 
treatment of the intersite coupling in a perturbation theory proved satisfactory to tackle a 
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lattice which nonhnearity is stronger than the dispersion. In contrast, when the interaction 
between first neighbors dominates the on-site nonhnearity, which remains however non-zero, 
the amphtude of the dynamical structure factor makes visible the biphonon, although the 
energy transfer is almost harmonic. In such a case, the variation of the biphonon DSF 
with respect to the transfer momentum q exhibits some peaks that are much larger than 
the multiple-scattering. This would hint that, provided that the DSF could be resolved 
accurately in g, the nonlinear behavior of certain materials could be worked out, even though 
their spectra show no apparent energy anharmonicity. 

In certain metal hydrides whose hydrogen anharmonicity is significant, we proposed a 
scheme to work out from the INS spectra, the potential landscape of the interstitial gaps. 
To approach some concrete cases, the theory requires, thought, further developments to con- 
sider, for instance, the three dimensionality of a realistic sample or the incoherent scattering. 
In a first step, this could be achieved within our perturbation theory. In addition to the 
possibility of simulating the incoherent INS in certain hydrides and molecular crystals, as 
those quoted in the present paper, it might be worth studying the nonlinear surface modes 
that could be investigated practically by different technics, e.g., electron scattering or in- 
frared adsorption. The low dimensionality of the surface could give an advantage to achieve 
a simulation that would be physically accurate, particularly upon the lattice geometry and 
various polarizations. As an example, we note the case of the CO molecules adsorbed on 
a Ru(lll) surface that has been studied both theoretically and experimentally^'^'^. As 
a low dimensional system, a vicinal surface, which might exhibit a regular one-dimensional 
nano-structure over several micrometers, would be the ideal substrate to explore the coherent 



scattering of a appropriate nonlinear surface modes. The behavior o 



a vicinal surface presents some specificity (see for instance Refs JsdiSTI and references therein) 



molecules adsorbed on 



and thus the feasibility of such a study remains under question. Though, the aim would 
be to carry out a momentum resolved experiment, as those attempted on the PtCl ethylene 
diamine chlorate^"^, to check whether the predictions established within the KG model 
are confirmed in practice, especially about the biphonon contribution to the DSF. Another 
possibility to test our theory, would be on the O2 solid^^i^ where the oxygen cross-section 
is mainly coherent. Although, in the /?-phase of the O2 crystal, the inter-molecular coupling 
overpasses the anharmonicity of the O2 stretch^^, our study showed that in such a case the 
inelastic scattering due to biphonon could yet be larger in magnitude than the two-phonon 
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response. Further, the phonon bound states could emerge at higher energy transfer as found 
for the triphonon and quadriphonon in metal hydrides^. At least from a theory viewpoint, 
the /3-phase of the O2 solid might be worth reexamining with recent experimental devices. 
Finally, according to the authors of Ref|3|, the potential landscape of hydrogen in PdH, that 
has been worked out from the INS spectral agrees well with a first-principles calculation 
made independently^. On the basis of that successful comparison, one may expect that 
a standard first-principles theory could be used to calibrate the KG parameters so as to 
simulate ab-initio the dynamical structure factor. 
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Figure captions 

Fig. [TJ Energy spectrum of a chain composed of unit cells, = 33 (full circles) and 
= 23 (empty circles): (a) the optical phonon branch for model parameters C = 0.05, 
^3 = 0.12 and ^4 = 0.01; (b) the two-phonon energy region for the same parameters as 
(a); (c) the same as in (b) but for C = 0.05, A3 = 0.05 and A4 = 0.01. The Y axis unit is 
hQ. The wave vector is reported on the X axis, whose unit is (vr/ao) and ranges over the 
lattice first Brillouin zone. 

Fig. [21 (Color online) A 3D plot of the inelastic S{q, u) as a function of the dimensionless 
energy transfer u and the scalar product q = [q.u] of the transfer momentum q and the 
polarization u. The parameters are the same as in Figs. [T] (a) and (b). The figure is easier 
to depict in color. 

Fig. [3l A plot of the Debye-Waller factor [S{q,uj = 0)/A^] versus the scalar product 
q = [q.u], for different model parameters. Each symbol corresponds to a Bragg peak. Our 
numerical results are reported for C = 0.05, A3 = A4 = (diamonds) and C = 0.05, A3 = 0, 

= 0.2 (triangles). The exact formula derived in a harmonic model (dashed line) confirms 
the former while the latter has been computed within a perturbation theory (circles) too. 
The X axis, which unit is Tr/ag bears the scattering vector q over a range of 60 Brillouin zones. 

Fig. m A plot of the coherent inelastic S{q,uj) profile versus q = [q.u], for phonon in 
a harmonic lattice: C = 0.05, ^43 = ^44 = 0. The thin dashed line has been obtained by 
a perturbation theory. The thin solid line corresponds to our numerical treatment and 
the thick dashed line to a plot of the exact formula^ (the two latest curves are hardly 
distinguished in the graph). In the inset, the area in the dot-dashed rectangle is magnified. 
The X axis unit is vr/ao. 

Fig. [5l A plot of the coherent inelastic S{q,u!) profile versus q = [q.u], for phonon, 
biphonon and two-phonon states. The model parameters are: (a) C = 0.05, ^3 = 0.12 and 
A4 = 0.01; (b) C = 0.05, A3 = 0.05 and A4 = 0.01. The dashed lines correspond to our 
perturbation theory and the solid lines to our numerical treatment. The X axis unit is vr/ao 
and the scattering vector q ranges over 60 Brillouin zones. The inset shows a magnification 
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of the two-phonon contributions. 

Fig. [HI In (a), the same as in Fig. [S](a) but for different model parameters: C = 0.06, 
^3 = 0.17 and A4 = 0.0231. In addition, have been plotted the S{q,uj) profile for the 
triphonon and the integral of S{q,uj) over the energy transfer u (dashed line). The elastic 
contribution has been skipped. In (b), a plot of the profile of the S{q,u!) integral over the 
reciprocal lattice (see the text), versus the energy transfer uj, for same parameters as in 
(a). The inset in (b) shows the corresponding potential landscape V{Xj) and its quantum 
levels. Our energy unit HQ has been fixed to 107 meV. 
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